#include <math.h>

/* flip() and clean() are supporting routines called by sieveOfAtkin()
 */
int _sieve_flip(int list[], int n) {
    if (list[n] == 1) list[n] = 0;
    else list[n] = 1;
    return 0;
}

int _sieve_clean(int list[], int n, int upto) {
    int k, l = n*n, i;
    for (i = 2, k = l; k <= upto; i++, k = i * l) list[k] = 0;
    return 0;
}

/* Generates prime numbers upto the integer specified by primes_upto.
 * Stores the result in the array primes[] and returns number of primes
 * generated upto primes_upto.
 */
int sieveOfAtkin(int primes[], int list[], int primes_upto) {
    int i, j, x, y, d_limit, n;

    list[0] = list[1] = 0 ;
    list[2] = list[3] = 1 ;
    
    d_limit = (int) (sqrt((float) primes_upto));

    for (i = 4; i <= primes_upto; i++) list[i] = 0;

    for (x = 1; x <= d_limit; x++) {
        for (y = 1; y <= d_limit; y++) {
            n = 4 * (x * x)+(y * y);
            if (n <= primes_upto && (n % 12 == 1 || n % 12 == 5)) _sieve_flip(list, n);

            n = 3 * (x * x)+(y * y);
            if (n <= primes_upto && n % 12 == 7) _sieve_flip(list, n);

            n = 3 * (x * x)-(y * y);
            if (x > y && n <= primes_upto && n % 12 == 11) _sieve_flip(list, n);
        }
    }

    for (n = 5; n <= d_limit; n++) {
        if (list[n]) _sieve_clean(list, n, primes_upto);
    }

    for (i = 2, j = 0; i <= primes_upto; i++) {
        if (list[i]) primes[j++] = i;
    }

    return j - 1;
}
